/*
 * matrix.cpp
 *
 *  Created on: 04/06/2011
 *      Author: Yotam and Gilad
 */

#include <iostream>
#include "mtm.h"
#include <cassert>
#include "matrix.h"

using namespace std;
using namespace MatrixExceptions;

/**
 * squareDimension: checks if a given matrix is a square. if it is, returns its
 * 	dimension. if it isn't, returns 0
 */
int squareDimension(const Matrix& matrix) {
	return (matrix.columns() == matrix.rows()) ? matrix.rows() : 0;
}

void Matrix::verifyIndex(int row, int col) const {
	if ((row <= 0) || (row > numberOfRows) || (col <= 0) || (col
			> numberOfColumns)) {
		throw IndexOutOfRange();
	}
}

void Matrix::compareDimensions(int a1, int a2, int b1, int b2) {
	if ((a1 != a2) || (b1 != b2)) {
		throw BadDimensions();
	}
}

double* Matrix::allocateAndCopy(const double* data, int size) {
	double* newData = new double[size];
	for (int i = 0; i < (size); ++i) {
		newData[i] = data[i];
	}
	return newData;
}

Matrix::Matrix(const int& rows, const int& cols) :
	numberOfRows(rows), numberOfColumns(cols) {
	if ((rows <= 0) || (cols <= 0)) {
		throw BadDimensions();
	}
	data = new double[rows * cols];
	for (int i = 0; i < (rows * cols); ++i) {
		data[i] = 0.0;
	}
}

Matrix::Matrix(const Matrix& matrix) :
	numberOfRows(matrix.numberOfRows), numberOfColumns(matrix.numberOfColumns),
			data(allocateAndCopy(matrix.data, numberOfRows * numberOfColumns)) {
}

Matrix::~Matrix() {
	delete[] data;
}

int Matrix::columns() const {
	return numberOfColumns;
}
int Matrix::rows() const {
	return numberOfRows;
}

const double& Matrix::operator()(const int& row, const int& col) const {
	this->verifyIndex(row, col);
	return data[(row - 1) * numberOfColumns + (col - 1)];
}

double& Matrix::operator()(const int& row, const int& col) {
	this->verifyIndex(row, col);
	return data[(row - 1) * numberOfColumns + (col - 1)];
}

bool Matrix::operator==(const Matrix& matrix) const {
	compareDimensions(numberOfRows, matrix.numberOfRows, numberOfColumns,
			matrix.numberOfColumns);
	for (int i = 0; i < (numberOfRows * numberOfColumns); ++i) {
		if (!(IS_DOUBLE_EQUAL(data[i], matrix.data[i]))) {
			return false;
		}
	}
	return true;
}

bool Matrix::operator!=(const Matrix& matrix) const {
	return !(this->operator==(matrix));
}

Matrix& Matrix::operator=(const Matrix& matrix) {
	if (this == &matrix) {
		return *this;
	}
	double* newData = allocateAndCopy(matrix.data, (matrix.numberOfRows
			* matrix.numberOfColumns));
	delete[] data;
	numberOfRows = matrix.numberOfRows;
	numberOfColumns = matrix.numberOfColumns;
	data = newData;
	return *this;
}

Matrix& Matrix::operator+=(const Matrix& matrix) {
	this->compareDimensions(numberOfRows, matrix.numberOfRows, numberOfColumns,
			matrix.numberOfColumns);
	int dimensions = numberOfRows * numberOfColumns;
	for (int i = 0; i < dimensions; ++i) {
		this->data[i] += matrix.data[i];
	}
	return *this;
}

Matrix& Matrix::operator-=(const Matrix& matrix) {
	return this->operator+=(-matrix);
}

Matrix& Matrix::operator*=(const Matrix& matrix) {
	this->compareDimensions(numberOfColumns, matrix.numberOfRows);
	int newRows = numberOfRows, newCols = matrix.numberOfColumns;
	int multipliers = numberOfColumns; // also = matrix.numberOfRows
	Matrix newMatrix(newRows, newCols);
	for (int i = 1; i <= newRows; ++i) {
		for (int j = 1; j <= newCols; ++j) {
			for (int k = 1; k <= multipliers; ++k) {
				newMatrix(i, j) += this->operator()(i, k) * matrix(k, j);
			}
		}
	}
	this->operator=(newMatrix);
	return *this;
}

Matrix& Matrix::operator*=(const double& scalar) {
	Matrix newMatrix = *this;
	for (int i = 1; i <= numberOfRows; ++i) {
		for (int j = 1; j <= numberOfColumns; ++j) {
			newMatrix(i, j) *= scalar;
		}
	}
	this->operator=(newMatrix);
	return *this;
}

Matrix Matrix::operator+(const Matrix& matrix) const {
	Matrix sumMatrix = *this;
	return sumMatrix += matrix;
}

Matrix Matrix::operator-(const Matrix& matrix) const {
	return this->operator +(-matrix);
}

Matrix Matrix::operator-() const {
	return (*this) * -1.0;
}

Matrix Matrix::operator*(const double& scalar) const {
	Matrix newMatrix = *this;
	return newMatrix *= scalar;
}

Matrix Matrix::operator*(const Matrix& matrix) const {
	Matrix newMatrix = *this;
	return newMatrix *= matrix;
}

Matrix operator*(const double& scalar, const Matrix& matrix) {
	return matrix * scalar;
}

Matrix Matrix::createUnit(const int dimension) {
	Matrix unitMatrix(dimension, dimension);
	for (int i = 1; i <= dimension; ++i) {
		unitMatrix(i, i) = 1.0;
	}
	return unitMatrix;
}

Matrix Matrix::transpose() {
	Matrix newMatrix(this->numberOfColumns, this->numberOfRows);
	for (int i = 1; i <= this->numberOfRows; ++i) {
		for (int j = 1; j <= this->numberOfColumns; ++j) {
			newMatrix(j, i) = this->operator ()(i, j);
		}
	}
	return newMatrix;
}

void Matrix::divideSubtractSwap(int row, int column, const double& scalar,
		Matrix& solutionMatrix) {
	devideRowByScalar(row, scalar); // make the leading element 1
	solutionMatrix.devideRowByScalar(row, scalar);
	for (int i = 1; i <= numberOfRows; ++i) {
		if (i != row) { //for each row but the subtracter
			double elementToZero = this->operator()(i, column);
			subtractRows(i, row, elementToZero);
			solutionMatrix.subtractRows(i, row, elementToZero);
		}
	}
	swapRows(row, column); // move the row to its right place
	solutionMatrix.swapRows(row, column);
}

void Matrix::swapRows(int row1, int row2) {
	for (int j = 1; j <= numberOfColumns; ++j) {
		double temp = this->operator ()(row1, j);
		this->operator()(row1, j) = this->operator ()(row2, j);
		this->operator()(row2, j) = temp;
	}
}
void Matrix::devideRowByScalar(int row, const double& scalar) {
	for (int j = 1; j <= numberOfColumns; ++j) {
		this->operator()(row, j) = this->operator()(row, j) / scalar;
	}
}
void Matrix::subtractRows(int row1, int row2, const double& scalar) {
	for (int j = 1; j <= numberOfColumns; ++j) {
		operator()(row1, j) = operator()(row1, j) - (operator()(row2, j)
				* scalar);
	}
}

Matrix Matrix::inverse() const {
	int dimension = squareDimension(*this);
	if (!dimension) {
		throw BadDimensions();
	}
	Matrix matrixCopy = *this;
	Matrix solutionMatrix = Matrix::createUnit(dimension);
	double leadingElement = 0.0;
	for (int column = 1; column <= dimension; ++column) {
		bool leadingFound = false;
		for (int row = column; row <= dimension; ++row) {
			if (!IS_DOUBLE_EQUAL(0.0, matrixCopy(row,column))) {
				leadingFound = true;
				leadingElement = matrixCopy(row, column);
				matrixCopy.divideSubtractSwap(row, column, leadingElement,
						solutionMatrix);
				break;
			}
		}
		if (!leadingFound) {
			throw Singular();
		}
	}
	return solutionMatrix;
}

